Fluctuating and dissipative dynamics of dark solitons in quasi-condensates 
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The fluctuating and dissipative dynamics of matter-wave dark solitons within harmonicaUy 
trapped, partially condensed Bose gases is studied both numerically and analytically. A study 
of the stochastic Gross-Pitaevskii equation, which correctly accounts for density and phase fluctua- 
tions at finite temperatures, reveals dark soliton decay times to be lognormally distributed at each 
temperature, thereby characterizing the previously predicted long lived soliton trajectories within 
each ensemble of numerical realizations (S.P. Cockburn et al, Phys. Rev. Lett. 104, 174101 (2010)). 
Expectation values for the average soliton lifetimes extracted from these distributions are found to 
agree well with both numerical and analytic predictions based upon the dissipative Gross-Pitaevskii 
model (with the same ab initio damping). Probing the regime for which O.SfesT < /x < l.GfesT, we 
find average soliton lifetimes to scale with temperature as r ~ T"'', in agreement with predictions 
previously made for the low-temperature regime fcsT ^ /i. The model is also shown to capture 
the experimentally-relevant decrease in the visibility of an oscillating soliton due to the presence of 
background fiuctuations. 

PACS numbers: 67.85.-d, 03.75.Lm, 05.45.Yv 



I. INTRODUCTION 

As intriguing realizations of quantum objects on a 
macroscopic scale, the existence of solitons within atomic 
Bose-Einstein condensates (BECs) has led to much ex- 
perimental and theoretical work. Both dark [IHH] and 
bright solitons [9l-fT2] have been observed within single 
species BECs. The existence of each can be argued as 
being linked to the robustness of the system's nonlinear 
waves and ultimately with the system's (near-) integra- 
bilty. This is chiefly responsible for the absence of dissi- 
pation, for example, during soliton-soliton collisions [5]. 
For a decay mechanism to manifest, a substantial depar- 
ture from the integrable limit must be enforced. 

BEC experiments typically take place within an effec- 
tively harmonic three-dimensional (3D) trapping poten- 
tial, the introduction of which already lifts the integra- 
bility of the system. In principle, this means that soli- 
tons are no longer protected from decay by the infinite 
number of conservation laws, as in the integrable homoge- 
neous ID Gross-Pitaevskii equation (GPE). In order that 
solitons be rendered dynamically stable against decay in- 
duced by transverse excitations [4j, they should be pro- 
duced in highly elongated gases [H] , in which phase fluc- 
tuations are known to play an enhanced role |14| . Then, 
in the absence of a thermal cloud, they are found to be 
stable in the special case of a longitudinal harmonic con- 
fining potential !15J , which has been shown to be a direct 
consequence of the periodic emission and reabsorbtion of 
sound waves as the solitons oscillate in a harmonic po- 
tential in [17]. 



rather rapidly upon reaching the condensate edge, an ef- 
fect attributed to thermal decay [131 dH]- More recent 
experiments have produced solitons with much longer 
lifetimes allowing for the observation of head-on colli- 
sions between solitons [6l [19] and the dynamics of one 
or more soliton oscillations [S] |B] . In Ref. [5 , dark soli- 
tons were created by a phase imprinting technique and 
found to exist for very long times, with a clear oscillatory 
trajectory of the density notch visible in the experimen- 
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The first successful observations of dark solitons in 
BECs were made by Burger et al. [1] in a somewhat elon- 
gated set-up, although the solitons were found to decay 
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FIG. 1. Illustrative dissipative dark soliton dynamics in the 
presence (upper plot) or absence (lower plot) of stochastically- 
sampled fluctuations for T « 1.9T,/, « O.STqc; Dissipative 
effects lead to the decay of the soliton in each case, while 
fluctuations (in the stochastic simulations) additionally lead 
to reduced soliton visibility and less regular oscillations. 



tal data. Interestingly, a strong shot-to-shot variation in 
soliton Ufetimes was also reported in this work: solitons 
were found to exist in single realizations for times an or- 
der of magnitude longer than an average trajectory was 
obtainable experimentally [3 [2D]. While attributed to 
small preparation errors in this experiment, phase fluc- 
tuations are typically important in quasi one-dimensional 
(ID) systems, suggesting that a similar effect might be 
expected to occur when introducing a soliton within a 
phase fluctuating condensate [211 122] ■ 

Given the competition between fluctuations and dissi- 
pation observed in experiments on dark solitons, in this 
paper we perform a detailed numerical and analytical 
study of the effect of each of these on the motion of dark 
solitons within harmonically trapped, phase fluctuating 
BECs. Figure l] shows an example of the stochastic (up- 
per panel) and dissipative (lower panel) dynamics which 
we observe via the theories we analyze here, discussed in 
Sec. II. 
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FIG. 2. (Color online) Form of the equilibrium dissipation 
7(2) at T = 250nK when the SGPE density (black solid) and 
the Thomas-Fermi density (blue dot-dashed) are used in the 
mean-field potential U{z). The vertical dotted lines indicate 
the T = Thomas-Fermi radius, R. 



II. THEORETICAL APPROACH 

Given the particle-like behavior already observed for 
dark solitons experimentally [S] [51 [IS]j the notion of a 
'heavy' soliton oscillating within a background of lighter 
thermal particles makes it tempting to draw an anal- 
ogy to the Brownian motion of a particle moving within 
a fluid of lighter particles, undergoing many scattering 
events. The interaction between a dark soliton and ther- 
mal excitations in a BEC was considered in [23, by means 
of a kinetic equation, which was found to be of Fokker- 
Planck form under the assumption that the momen- 
tum transfer per soliton-excitation interaction is much 
smaller than the soliton momentum. In this work, we 
adopt a complementary Langevin approach which should 
thus also be well suited to the study of such systems, 
having been originally conceived for just this purpose 
[231 [53], and we now discuss the pertinent model for 
weakly-interacting gases, known as the stochastic Gross- 
Pitaevskii equation (SGPE). 



A. Stochastic Gross-Pitaevskii equation 

Since highly elongated trapping potentials are neces- 
sary if a soliton is to be stable against transverse insta- 
bilities, the inclusion of phase and density fluctuations 
is likely to prove essential in capturing all the salient 
aspects of such necessarily "low dimensional" systems. 
The SGPE [251 - f25] . which we choose to employ here, is 
a model well suited to this task as it captures both the 
dissipative and fluctuating dynamics inherent to finite 
temperature BECs, while additionally satisfying the re- 
quired balance between these two factors. 

Within this scheme, the low-energy modes of our effec- 
tively ID system may be represented by the stochastic 



differential equation 
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^{z,t)+rjiz,t), 



+ V{z) 



(1) 



where ^(2;,t) is a complex parameter describing the oc- 
cupation of such low-lying modes, V{z) = {l/2)muj'^z'^ is 
the axial trapping potential, g — 2huj±a is the effective 
ID interaction strength [29 (with a the s-wave scatter- 
ing length), and 77(2, t) is a complex Gaussian noise term, 
with correlations given by the relation {rj* (z, t)ri{z' , t')) = 
2/17(2, t)kBT6{z - z')S{t - 1'). The strength of the noise, 
and damping, due to contact with higher energy thermal 
modes, is given by 7(2:, t). 

While the soliton will reside in the low-lying modes, 
in arriving at Eq. ^ it is assumed that the high energy 
thermal atoms in the system may be treated as though at 
equilibrium; this should be valid for small perturbations, 
e.g. when introducing a dark soliton into a large system 
{R ^ £,, where ^ = h/y/frmg essentially sets the soliton 
width) which would have little effect on the high-lying 
modes, which should instead remain close to equilibrium. 
In this picture, the high- and low-energy systems are as- 
sumed to be in diffusive and thermal equilibrium with a 
common temperature T and chemical potential /i. More- 
over, due to the large mode occupations typical within 
degenerate Bose gases, we are well justified in the fur- 
ther assumption that the low-energy modes are highly 
occupied for the classical field approximation to be valid 



B. Form of Dissipation 

Within the formulation of Stoof [5S1 [17] , the quantity 
which parametrizes the strength of the noise and damp- 
ing is the so-caUed Keldysh self-energy /il]^(z,i), which 
may be related to the dissipation term in Eq. (fTj) via 
7(z,t) = il3hT,^{z,t)/4: [171I3T]- Following the methods 
of Refs. [211 [33, and assuming the system to be close to 
equilibrium, the integral to be evaluated for the dissipa- 
tion is 

l{z) = ^^^_^ / rf^2 / dEa S{h,k2,k3){l + Ni)N2N3 

(2) 

where < E^ < 00, N, = [exp{l3{Ei + U{z)- fi))-l]-^ is 
a Bose-Einstein distribution, U{z) = V{z) + 2g{\'^{z)\'^), 
El = E2 + E3 + U{z) - ^, and 

S{ki,k2, ks) =- {sgn(fci + k2 - k^) + sgn(fci - fc2 + fcs) 

-sgn(A:i + /S2 + k^) - sgn(fci - ^2 - ^3)} , 

(3) 

with ki = ^J^/mEijW . 

The damping may be evaluated numerically (e.g. us- 
ing Gauss-Legendre quadrature), an example of which 
is shown in Fig. [2j The data shown is for T = 250nK 
and illustrates two cases, when the density in the poten- 
tial U{z) is given by the SGPE density (black solid line) 
or the T = Thomas-Fermi density (blue dot-dashed 
line) [cf. Fig. 2 of Ref. [3T]]. A notable feature in the 
form of 7(2;) is the peak close to the Thomas-Fermi ra- 
dius, which indicates that the scattering rate is largest 
near the edge of the quasi-condensate. Also, the peaks 
in 7(2:) calculated using the SGPE density (black solid 
curve) are noticeably closer to the trap center than for 
the T = Thomas-Fermi case at the same [i (blue dot- 
dashed curve), due to the effects of thermal depletion on 
the equilibrium density profile. We also note that the 
shape of 7(2) is qualitatively similar to the spatial form 
of the scattering rates found within numerical implemen- 
tations of the Zaremba-Nikuni-Griffin (ZNG) approach 

ElEZ]. 



C. Reduction to a dissipative Gross-Pitaevskii 
equation 

Neglecting the noise term of Ec^. (IT]) leads to the fol- 
lowing dissipative GPE (DGPE) for the condensate wave- 
function (see e.g. the review of Ref. 38 ): 



,;,MM = (i_,^(,,i)) 



dt 



2m dz^ 



Viz) 



This is expected to be a good model when damping ef- 
fects are dominant over diffusion, and representative of 
the mean field soliton dynamics in this case. An advan- 
tage of deriving this model from the SGPE is that the 
damping may be calculated in an ab initio manner, via 
Eq. (l2|, rather than selected phenomenologically [39ll40j . 



D. Overview of numerical procedure 

Before studying the soliton dynamics, an important 
preliminary step is to generate a suitable equilibrium 
state. This state differs between the DGPE and SGPE 
approaches: for the SGPE, the initial condition (the 
noisy field denoted by ^ in Eq. (Ilj) is dynamically ob- 
tained at each temperature by equilibration with the 
higher energy heat bath atoms - see, e.g., Ref. [571 EH" 
1331 . For the DGPE, we instead use the ground state of 
the GPE (the complex field denoted by (j>), which is ob- 
tained efficiently by imaginary time propagation of the 
GPE. We make the latter choice because the imaginary 
time solution to the GPE is also the ground state of the 
DGPE, the action of which is to drive the system to- 
wards this state, at a rate related to 7. If instead we were 
to start with the SGPE initial condition as an input to 
the DGPE, this would mean the system was already out 
of equilibrium, even prior to introducing a perturbation 
such as a dark soliton. 

To generate a dark soliton in both stochastic and purely 
dissipative simulations, we multiply the equilibrium state 
by the soliton wavefunction 



Aoiiz) = (tanhiCz/S,) + i{v/c), 



(5) 



+ ffl<^l'-A' 



(4) 



(/.(Z,t). 



where ( — -^/l — (w/c)^, ^ is the healing length, and c is 
the speed of sound. The input soliton velocity, Uinput, is 
chosen here such that |uinput| = 0.25c in order to gener- 
ate a relatively deep density notch (which is thus clearly 
identifiable over the thermal background noise). We fo- 
cus on this method of initial state preparation, rather 
than full phase imprinting, as we wish to highlight the 
role of both the initial and dynamical thermal noise even 
for a completely reproducible preparation mechanism. 

As (f>{z) of Eq. Q is a mean field, only one realization 
of the DGPE is required at each temperature. For the 
SGPE simulations, we must instead generate an ensemble 
of initial states, {^}, consisting of several hundred real- 
izations, which we initially propagate to equilibrium with 
the higher energy thermal cloud. The equilibrium state 
is parametrized by the chemical potential, fi, and tem- 
perature, T, and both fi and 7(2; /x,T) remain the same 
between the SGPE and DGPE simulations (the temper- 
ature dependence in the DGPE arises only through 7, 
while the strength of the noise in the SGPE is also tem- 
perature dependent). 

We choose to work with realistic trapping frequencies 
ujz = 27r X 10 Hz and u}j_ — 2t: x 2500 Hz, and set 
ji = 395/1^2, which corresponds to around 20000 *^Rb 
atoms (at T — Q). To strike a balance between practical 



simulations and interesting dynamics due to the interplay 
between dissipation and fluctuations, we focus here on a 
temperature range T = 150-300nK, which corresponds 
to 0.16 < T/Tqc < 0.34 and 6 < T/T^ < 13, where 
Tqc = N{nuj^)y\n{2N)kB M and T^ = 7V(fej^)VMB 
[l4| respectively define 'characteristic' temperatures for 
the onset of density and phase fluctuations. Hence we 
deal with a partially condensed Bose gas, which is well 
within the phase fluctuating regime. 



III. ANALYSIS OF STOCHASTIC DYNAMICS 

The link between the quantization of classical theo- 
ries permitting soliton solutions, and dissipative quan- 
tum systems was discussed in Refs. I 15H47J . There, it was 
highlighted that just as the motion of a classical parti- 
cle within a viscous environment has both a damped and 
fluctuating component, the situation is the same in the 
quantum case. The motion of such a classical particle 
can then be characterized by two system properties, a 
damping due to the systematic force applied to the par- 
ticle and a diffusion related to this interaction [JS] . This 
was shown to be true also in the quantum case f49l , in 
which the dissipation manifests instead as a damping of 
the particle wavepacket center of motion, whereas diffu- 
sive effects lead to a spreading of the wave function for 
the particle [17]. For a soliton, the former would lead to 
a damping of the motion of the soliton center, while the 
latter would give rise to an increased uncertainty in the 
soliton position. 

For soliton solutions to integrable, classical one- 
dimensional theories, in which case dissipation has no 
role, the propagation in space is undamped and the soli- 
ton dynamics is essentially captured by knowledge of the 
soliton center. However, in the quantized field-theory at 
finite temperature, not all degrees of freedom "collabo- 
rate" in the formation of the soliton Wf] and the result is 
a residual interaction which is shown to lead to a Brown- 
ian type motion of the soliton. Therefore, the dynamics 
is no longer entirely captured by knowledge of the cen- 
ter of mass alone, as the diffusive nature also becomes 
important. 

In Ref. [17], the Brownian nature of the soliton mo- 
tion is related to the coupling of the soliton to the other 
system modes and excitations due to the presence of the 
soliton. Damping and diffusion in the quantum dissipa- 
tive system must be temperature dependent, since the 
excitations which scatter from the soliton are thermally 
activated. Returning to the BEC context, there is an 
obvious analogy between this work and a soliton propa- 
gating within a finite-temperature BEC. 

Studies of weakly-interacting, ID homogeneous Bose 
gases, considering the purely dissipative dynamics, found 
solitons to decay with a lifetime which varies with tem- 
perature as r - T""* for ksT < /x [23l[50]. For ksT » ^, 
it was predicted that r ^ T~^. (This temperature de- 
pendence was also found in a study on polarons [46) , and 



in the more general case of the "quantum impurity prob- 
lem" , applied in the setting of a heavy particle within 
a Luttinger liquid |51|.) The parameters in the present 
study are chosen such that 0.8/i < ksT < 1.6/i, meaning 
we probe an intermediate temperature regime which is 
hard to treat analytically yet in which soliton decay can 
occur on a convenient time scale. Such a regime has in 
fact been reached in a number of experiments on an atom 
chip, see e.g. Ref. [52] and references therein, although 
this type of set-up has yet to be used to observe the in- 
teresting soliton dynamics anticipated in such systems. 
Despite this, several recent experiments have nonethe- 
less provided motivation for modeling beyond purely de- 
terministic dissipative dynamics, with solitons observed 
to exist for times much longer than a reproducible av- 
erage trajectory could be produced 5 ; indeed prelimi- 
nary work on thermal decay with O.lfi < ksT < 0.5/i 
already found a significant spread in single trajectories 
j20j . and such an effect is expected to be amplified as 
fiuctuations become more important. The necessity for 
repeated runs in the SGPE formalism also has a strong 
link to the experimental approach of repeatedly creating 
successive BECs, so incorporates this effect naturally. 

Finally, we note that Ref. [50] discussed the fact that 
a full treatment of the soliton dynamics should satisfy a 
fiuctuation-dissipation theorem, which is not the case in 
retaining only the damping aspects of the thermal back- 
ground. Fluctuations however complicate the analysis of 
experiments and the stochastic simulations we have un- 
dertaken, as discussed in the following Section. 



A. Extracting soliton information 

Within the SGPE approach, observables are obtained 
by sufficient averaging over many stochastic realizations. 
For example, the density for an average over N noise re- 
ahzations, '^i,'i'2, ■ ■ ■ , ^jv = {^jw is given by {n{z)) = 
(^*(2)*(z)) = (l/iV)^,^,**(^)*,(z). In order to 
extract information on the soliton present within each 
stochastic realization, we find we must extract the soli- 
ton trajectory prior to performing any averaging. If we 
instead simply average over the set of density profiles to 
obtain (n(z)), the soliton is quickly found to be "washed 
out" , despite the fact that a soliton is still present within 
each individual realization ^Pi, ^2, ■ • ■ , ^jv ^ this may be 
understood from the fact that after some time solitons 
within different realizations are likely to be at different 
positions, so averaging over many such density profiles 
quickly leads to a loss of information on the individual 
soliton positions. 

This effect is illustrated for the case of an initially 
static soliton in Fig. [3] The average density is plotted 
in Fig.jsja), showing the 'average soliton' filhng up over 
time; the density at the center of the sample, where the 
soliton is initially positioned, is shown in Fig. P][b) and 
grows towards the equilibrium value, with a growth in 
time which goes as ex (1 — exp(— Ft)), with r(> 0) the 



growth rate. The experimental analogue of this effect 
was discussed in Ref. [20], where it was noted that it is 
indeed single-shot soliton runs that should be analyzed to 
measure decay, rather than averaged images, as the soli- 
ton contrast is quickly found to be smeared out in the 
latter. In addition, Dziarmaga analyzed the diffusion of 
a soliton due to quantum fluctuations |53], while here it 
is thermal fluctuations which cause the diffusion process 
(see also [UlIiaiSlHSZ])- 

While soliton depth, or contrast, is quickly lost in av- 
erage images due to statistical effects on the motion, we 
now look at the fate of the solitons residing within each 
single SGPE realization. This is illustrated in Fig. |4J 
where we instead consider a moving soliton. The soliton 
depth, risoi, is measured relative to the background den- 
sity, meaning, for example, rigoi equal to the background 
density would correspond to the static, black soliton of 
Fig. [3j As in the case of a static soliton, the moving soli- 
ton is quickly lost in the average density profile, as shown 
by the density snapshots of Figs. |4[a)-(b). However, if 
we look instead at Fig. |4]jc), then we see that each sin- 
gle stochastic realization in fact still contains a very deep 
soliton: this plot shows the result of measuring rigoi, the 
depth of the soliton relative to the background density, 
within each individual realization of the ensemble {^Iat. 
Fig. [ijjc) shows the average of this set of depths, {usoi)- 
This average depth is quite different to that obtained by 
averaging over the individual density profiles (i.e. calcu- 
lating {^*{z)'9{z)) as in Figs. [4](a)-(b)), as importantly 




the soliton is first located within a single realization (and 
this location generally varies from realization to realiza- 
tion) before information on its depth is extracted. Sim- 
ply averaging the density profiles of individual runs after 
some time leads to a smooth profile, as the soliton den- 
sity minimum in a particular run becomes outweighed by 
the higher number of runs in which there is no soliton at 
that particular spatial point. In our analysis, we there- 
fore follow the density notch associated with the soliton 
present within each realization, as found also to be exper- 
imentally most consistent when considering soliton decay 

We can only track a soliton until the point where it 
becomes indistinguishable over the background density 
fiuctuations which defines the soliton decay time (in an 
analogous manner to experimental observations); carry- 
ing this out for each individual realization allows us to 
obtain an ensemble of soliton trajectories and so a distri- 
bution of soliton dynamics {^Jat which we can then an- 
alyze. This means of extracting data from the stochastic 
simulations appears consistent with the idea of the soli- 
ton center as a good quantum dynamical variable /i6[ . 



B. Spread in initial soliton velocity 



As discussed in Section II D[ we introduce a soliton of 
prescribed velocity (|uinput| — 0.25c) in an identical man- 
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FIG. 3. (Color online) (a) Average density profile from an 
ensemble of stochastic simulations, with each realization con- 
taining an initially static, black soliton. Thermal diffusion 
of individual solitons leads to the initially very deep notch 
average density becoming gradually filled. Snapshots are 
shown at several times between the initial state (black solid 
line) to the final state (dot-dashed green line), (b) Growth 
of the average central density tends towards the equilibrium 
value, as the soliton is lost in the average density profile, even 
though the solitons are still present within individual runs 
(see Fig. |4|. 



FIG. 4. Typo in text in (c) within the figure 'soltion'. (Color 
online) (a) Average density profile containing a moving soliton 
(jfinputl = 0.25c); (b) as in (a), but focused on the soliton 
region, (c) Average of the soliton depths extracted from each 
individual stochastic realization - this remains nearly constant 
meaning a deep soliton remains within each individual run. 
Information on the soliton depths is lost in the average density 
profile due to the soliton density notches of each run being, in 
general, at different spatial points. This implies that soliton 
information should be extracted from each realization, prior 
to any statistical analysis. 
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FIG. 5. (Color online) Comparison of the average measured 
soliton velocity (black circles) to the input velocity (dashed 
brown line) for various temperatures. The error bars indicate 
the standard deviation of the measured velocities. 



ner within each simulation. Despite this, one might ex- 
pect a spread in initial velocities due to the the fluctuat- 
ing background density, n — |^P, and the formula relat- 
ing this to the soltion velocity, v/c = ^\ — risoi/ri. Fig.p^ 
shows the mean initial velocities from a few hundred re- 
alizations versus temperature, with error bars showing 
the standard deviation. This plot shows a tendency for 
the soliton velocity to be lower than the prescribed value 
(shown by the dashed horizontal line) for higher temper- 
atures, and consequently larger fluctuations. 



C. Influence of initial vs. dynamical noise 

In addition to a variation in initial velocities, fluctu- 
ations due to the dynamical noise term of the SGPE 
also affect the solitons throughout their lifetimes (see e.g. 
Ref. [SS] for related work on fluctuations during vortex 
decay). We use the decay time, i.e. the time at which 
our tracking algorithm can no longer distinguish a soliton 
over the background noise, as an observable to measure 
the soliton dynamics. 

To study the relative importance of initial versus dy- 
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FIC. 6. (Color online) Scatter plots showing decay times, 
scaled to the ensemble average (y-axis), vs. initial soliton 
velocity, scaled to the input velocity (x-axis); T — (a) 150 uK, 
(b) 175 nK, (c) 200 nK (d) 250 uK and (e) 300 nK. (f) shows 
the standard deviation in the decay times (solid blue line) and 
initial velocities (dashed green line) vs. T. 



namical noise, we construct a scatter plot (see Fig.l6[a)-(e 
)) showing, on the x-axis, the initial velocity, v(t = 0), 
scaled to the input velocity used for deterministic soli- 
ton creation, winput- On the y-axis, we show the decay 
time, T, scaled to the ensemble average (r). The ver- 
tical dashed red line indicates the velocity Winput, while 
the horizontal dashed line indicates the ensemble aver- 
age decay time, (r). At all but the lowest temperature of 
Fig.|6][a), there are more points to the left of the expected 
velocity (vertical dashed line) than the right, suggesting 
that the soliton generated within many of the realizations 
had a velocity which was lower than the input velocity 
of Eq. ^ (consistent with Fig. 5). 

Looking instead at the vertical trend, we see that at 
each temperature there are several solitons which exist 
for very long times relative to the ensemble average (hor- 
izontal dashed line). Two examples of such long-lived 
trajectories are shown in the lower, rightmost plots of 
Fig.[Zl 

The data of Fig.[6]also allows us to indicate whether it 
is the variation in initial velocities or the dynamical noise 
which has greatest influence on the variation in soliton 
dynamics. The relative variation due to each of these ef- 
fects is quantified in Fig. loFf) which shows the standard 
deviation for each quantity, scaled to the relevant mean. 
It is clear that the relative spread in the decay times is 
far greater than the spread in the initial velocities, indi- 
cating that the dynamical noise has a larger cumulative 
influence on the observed variation in the soliton dynam- 
ics. 

Fig. [7] summarizes the large variation in soliton dy- 
namics that is observed, despite the identical means of 
soliton generation employed. The top panel shows a his- 
togram of soliton decay times, with a selection of rep- 
resentative soliton trajectories in the panels below. The 
colour/position of the trajectories in the lower plots cor- 
respond to the highlighted bins of the histogram. 

In the case of the leftmost and middle plots, we see the 
trajectories are no longer oscillatory after some time, and 
instead become noisy. This is the point at which the soli- 
ton is lost, and defines the decay time (after which the 
graph simply represents numerical noise). The middle 
plots show trajectories for solitons from the bin contain- 
ing the mean decay time (5th bin from left, highlighted 
blue), which are typically very close to the trajectory ob- 
tained from the DGPE simulation at this temperature. 
Finally, in the rightmost plots, we consider a bin corre- 
sponding to decay times which are longer than the en- 
semble mean; their amplitudes are far below that of the 
purely dissipative trajectory, shown by the dashed green 
line, indicating that the noise can act in some cases to 
stabilize a soliton against decay, relative to the purely 
dissipative evolution at the same temperature. 







FIG. 7. (Color online) Histogram of stochastic soliton decay times (upper plot) and two representative sets of trajectories from 
the highlighted bins (lower plots). The solid lines in the lower plots show stochastic trajectories while the dashed green lines 
is the purely dissipative result. The left/middle/right most plots have decay times within the corresponding left/middle/right 
highlighted histogram bins. The trajectories from the middle highlighted bin have a decay time close to the mean value, and 
we find the trajectories to be qualitatively similar to that of the dissipative case. 



D. Distribution of soliton decay times 

An example of the distribution of decay times at T = 
200 nK is shown in Fig.[8J It is clear that the distribution 
of sohton decay times is non-Gaussian, and we instead 
find it to be well fitted by a lognormal distribution, 



P{t) 



T^Jl-K 



cxp 



(In T — m) 
2(t2 



(6) 



where m is the mean and a the standard deviation of In t. 
A lognormal model is often applied to a system in which 
decay is caused by random events, which causes decay 
at a rate proportional to the amount already present, 
so in other words a runaway process. For example, Kol- 
mogorov suggested that for anything which decays in this 
multiplicative way, the time to failure should follow a log- 
normal distribution [59]. 

This distribution has a long tail at large decay times, 
which reflects the finite probability of some long lived 
solitons within a given ensemble of simulations. Inter- 
estingly, extreme cases of long-lived solitons have been 
observed experimentally within single runs, while the 
typical, reproducible decay time for the experiment was 
around an order of magnitude lower [5]. That soliton 
decay times follow a lognormal distribution is a possible 
reason for this observed behavior. 

It is also interesting to note that a lognormally dis- 
tributed soliton amplitude, which we denote 5'(i;7,/i), 
would solve a stochastic equation of the form 

dS(t; 7, /i) — mS{t; 7, ^)dt + aS{t; 7, /i)(iVF, (7) 
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FIG. 8. (Color online) Histogram of soliton decay times for 
an ensemble of few hundred SGPE runs (solid black circles) 
and a fit to a lognormal distribution (dotted blue). 



which is the stochastic differential equation defining a ge- 
ometric Brownian motion (here, m is the average growth 
rate, a is the so-called volatility and dW denotes a Weiner 
process). As the soliton depth determines how far from 
the trap center the turning point of the motion is, it 
is therefore directly related to the amplitude of the os- 
cillations. In turn, as the decay time is determined by 
reaching a certain depth, then we expect the distribution 
in this variable to display a similar lognormal distribu- 
tion. This would imply that the amplitude of the soliton 
oscillations undergoes a lognormal random walk. 

Eq. ([7|) should also be compared to the subcritical soli- 
ton equation of motion derived in Ref. [21J , and also that 
discussed in Section [Tv] [Eq.(24)]. The latter two pre- 
dict soliton oscillations with an exponentially increasing 
amplitude envelope dependent upon the damping, and 



the solution to Eq. ^ would also correspond to such an 
exponential function in the limit that cr — >■ (or, equiv- 
alently, fluctuations are neglected). Conversely, fluctua- 
tions might be introduced retrospectively to the dissipa- 
tive model by allowing the damping to become a stochas- 
tic variable: a similar idea has been applied in describing 
the collisions between optical solitons fSU], where it was 
found that collisions could be described by a nonlinear 
Schrodinger equation perturbed by stochastic parameters 
obeying strongly non-Gaussian statistics. Interestingly, 
in this case the soliton amplitude was also found to be 
lognormally distributed. 



E. Decay times vs. T: Analyzing the distributions 

In order to extract meaningful quantities from our re- 
sults, we must analyze the distributions of decay times 
which we obtain from the stochastic simulations. To do 
so, we proceed by fltting a normalized histogram of decay 
times at each temperature to Eq. m. At each tempera- 
ture, we flnd that the fit matches the numerical data well, 
suggesting that the underlying mechanism of the soliton 
decay is multiplicative in nature. The results of fitting 
the decay times for all temperatures considered are shown 
in Fig. [9j the inset shows the same data with the x-axis 
on a log-scale. We can see from the main plot of Fig. [9] 
that the distributions of decay times become increasingly 
shifted towards the origin with increasing temperature, 
as dissipative effects reduce the soliton lifetimes. The in- 
set highlights that it is the logarithm of the decay times 
which is Gaussian distributed. From the inset it is clear 
also that the variance of ln(T) increases as temperature 
is increased. 

An analysis of the behavior of the decay time distri- 
butions with temperature is given in Fig. 10 To extract 



tion of Eq. ^, is given by (r) = exp [^ -I- cr^/2]. The 



variance is instead calculated as Var[T] 
1) cxp [2fi - 



(exp[aV2] 



the average soliton behavior, we consider first the ex- 
pectation value, (r) , which for the lognormal distribu- 




FIG. 9. (Color online) Lognormal probability distributions 
obtained from fitting soliton decay time histograms obtained 
at the temperatures indicated. Inset: Same data with a log- 
scale X-axis; it is the logarithm of the decay time which is 
Gaussian distributed. 



cr^], and we use y^VarjrJ to generate the er- 
ror bars in the SGPE data of Fig. [lOja). 

Referring to Fig. [Ma), we flnd that (r) for the SGPE 
results (black circlesjclosely follows a T~^ behavior [55] , 
shown by the brown dashed line. Therefore the average 
behavior follows the same scaling predicted by the models 
in [531 EOj for fcsT <C ^. We work instead in the regime 
where fc^T ~ /i, which is difficult to treat analytically, 
however find the low temperature result to extend to this 
regime as well (although our numerics also relies on the 
classical field approximation). 

In order to compare meaningfully between the cur- 
rent SGPE analysis and the DGPE soliton dynamics, 
we must account for the changing level of background 
density fluctuations as temperature is varied within the 
SGPE [51] [5S]. To do so we have measured the minimum 
depth at which solitons can still be resolved in the SGPE 
runs, at each temperature [ST]. This is then used to ex- 
tract a soliton decay time within the DGPE simulations, 
i.e. the time for the soliton to decay to the tempera- 
ture dependent depth extracted from the SGPE simula- 
tions. The numerical DGPE results are shown by the 
fllled red diamonds in Fig. |10[ a), and display very good 
agreement with the ensemble average stochastic results, 
(t). We additionally show the results of our analytic 
model for the DGPE dynamics (hollow green squares), 
which is discussed in the following Section IV The de- 
cay times obtained within this model also agree well with 
the numerical DGPE and SGPE results, when limits on 
the soliton visibility due to thermal fluctuations are ac- 
counted for (see Section Iv] for details). 

As indicated above, the SGPE decay time distribu- 
tions are not symmetric (see Fig. [o]) and feature a long 
tail at long soliton lifetimes for all temperatures consid- 
ered. Our results suggest also that the noise stabilizes a 
certain number of solitons created within each stochas- 
tic ensemble of realizations against decay, relative to a 
soliton undergoing purely dissipative dynamics under the 
DGPE. This was also apparent from the decay time his- 
tograms shown in Refs. [5T1[SS]. The inset to Fig. floFa) 
shows the cumulative distribution function, D(t), which 
measures the fraction of solitons which have decayed at 
any time, based on the P{t) curves of Fig. [o] Clearly, 
dissipative effects are more dominant in the high tem- 
perature data, with all solitons found to have decayed at 
relatively short times. That the cumulative distribution 
function has a slow asymptote towards unity, is another 
way to see that some fraction of the solitons live for far 
longer than the average decay time. It is also interesting 
to note that scaling the time axis in each curve to the 
average at that temperature (r), we flnd each graph to 
collapse down close to a single curve. 



Fig. 10 ^b) compares several different measures of the 
distribution function: the arithmetic mean, or (r), the 
geometric mean and the modal value. Comparing these, 
we see that the modal values give consistently lower decay 
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FIG. 10. (Color online) Analysis of the decay time distributions: (a) Expectation value of the SGPE soliton decay time (black 
circles; error bars denote one standard deviation) versus temperature. Corresponding DGPE decay times are also shown (filled 
red diamonds) and a fit to a function ex T~* (brown dashed line). Green squares show analytical predictions from the analysis 
of Section [rV] the inset shows the cumulative probabilty distribution: the fraction of solitons which have decayed, versus time, 
for each temperature, (b) SGPE expectation value (solid black) compared to the geometric mean (brown solid line) and modal 
(green dot-dashed line) values for each temperature; (c) standard deviation and (d) skewness of r versus temperature. 



times than the expectation value, which is a consequence 
of the peak in the distributions of Fig. |9] being located to 
the left of the mean. 



In Fig. 10 'c) we measure the standard deviation of t, 
scaled to (r), which shows that the relative spread of 
the distributions increases monotonically with T. This 
is as one might expect, given that higher temperatures 
give rise to more fluctuations and so a wider variation 
between realizations. This also shows that while the low- 
est temperature distribution of Fig. [9] appears far wider 
than the highest temperature case, the scaled standard 
deviation (s.d.) illustrates that it is actually smaller. Fi- 
nally, Fig. 10 'd) shows the skewness of the distribution 



functions, which is also found to increase with tempera- 
ture. 

Next, we choose to decouple the effects of fluctuations 
and damping, and consider the dissipative behaviour 
alone through the mean field DGPE. A crucial ingredient 
in both the analytic and numerical calculations which we 
present, is the form of the damping term which we obtain 
from the SGPE formalism. 



which is the sum of two Lorentzians, centred near the 
edge of the quasi-condensate. The parameters a, c, and 
d depend in general on temperature, T, the trapping 
frequencies, chemical potential, and atomic species (al- 
though all except T remain fixed in our analysis) . Various 
values of these parameters, indicating their temperature 
dependence, are provided in Table |T] 

We will make use of the form given by Eq. (Is]) in the fol- 
lowing Section, in which we generalize our previous work 
on spatially constant dissipation reported in Ref. |21) , by 
developing analytic solutions for the more general case of 
spatially varying dissipation. 



A. Analytical approximations 

To perform our analytical work based on perturbation 
theory for dark solitons we follow the "routine" proce- 
dure of rescaling the equation in appropriate units. In 
particular, the ID dissipative GPE (DGPE) model may 



IV. ANALYTICAL APPROACH FOR 
DISSIPATIVE DYNAMICS 

We now consider the DGPE of Eq. Q. In order to 
proceed with analytical calculations, we first character- 
ize the form of 7(2:) in terms a simple analytical fitting 
function. We find the dissipation 7(2;) to be well approx- 
imated at all temperatures considered by the function 



7(z) 



(z + c)2+rf (z-c)2+(i' 



(8) 



T(nK) 


a 


c 


d 


125 


0.0137 


26.86 


2.555 


150 


0.0219 


26.63 


3.128 


175 


0.0330 


26.43 


3.464 


200 


0.0511 


26.29 


4.838 


250 


0.1031 


25.93 


7.114 


300 


0.1812 


25.69 


9.860 



TABLE I. Values of the parameters a, c and d [cf. Eq. 
for various values of temperature T. 



10 



0.008 




In the absence of the perturbation [P{v) — 0] , Eq. ( 10 1 



FIG. 11. (Color online) Shown is a comparison of 7(2;) as given 
by Eq. ^ [solid black line], the fit given by Eq. Q [dashed 
brown line], and the approximation of Eq. ( |15[ ) [dot-dashed 
blue line] for T — 150 nK. The dotted (red) curve depicts 
the mean value of 7(2;), while the vertical dotted green lines 
indicate the points z = ±i?, which define the rims of the TF 
cloud. The inset shows a magnification of this figure for the 
spatial interval [—R/2,R/2]: there, the agreement between 
the fit to 7(2:) (Eq. (|8|; dashed brown line) and the approxi- 
mate form of this (Eq. ( 15 1; dot-dashed blue line) is excellent. 



as the corresponding curves essentially coincide and become 
indistinguishable. 



is a conventional defocusing NLS equation, which pos- 
sesses a dark soliton solution of the form 1 62 : 



V 



'{z,t) — cos if tanh Z + i sin if, 



(12) 



where Z = cos (p[x — (sin (p)t\, and if is the "soliton phase 
angle" {\(p\ < 7r/2) describing the darkness of the soliton 
through the expression [ul^ = 1 — cos^ ipaech Z; note that 
the limiting cases ip = and cos t/s <C I correspond to the 
so-called black and grey solitons, respectively [S31 [M] . 
The effect of perturbation (111 on the dark soliton 



will be treated analytically by means of the adiabatic 
approximation of the Hamiltonian approach of the per- 
turbation theory for dark solitons. This approach was 
introduced in Ref. ^5\ for the case of a constant back- 
ground density, and subsequently used for trapped BECs 
in Ref. [66] (see also the review of Ref. [64] and our pre- 
vious work [H]). According to this approach, the pa- 
rameters of the dark soliton ( 12 ) become slowly- varying 
functions of time t, but the soliton's functional form re- 
mains unchanged. Thus, the soliton phase angle becomes 
if — > ip{t) and, as a result, the soliton coordinate becomes 
Z ^ Z = cos (p{t) [z — zo{t)], where 



also be written in the following dimensionless form, 
'1 



[i - -fiz)]dtTl' 



-di + v{z) + \^j\' - fi 



V-, (9) 



where the density jV'P) length, time and energy are re- 
spectively measured in units of 2a, a± = ^/H/mu>±, wj 
and huj±. In the case of a harmonic trap, the exter- 
nal potential takes the form V{z) — (1/2)J7^2:^, where 
f2 = u}z/u!± ^ 1 is the normalized trap strength, which 
is a naturally occurring small parameter of the system. 
The function 7(2:), which accounts for the dissipation, 
takes the form given in Eq.(|8|. 

We now seek a solution of Eq. ^ in the form -0(z, t) = 
'tjjb{z,t) eKp[—i9{t)]v{z,t), where ^i,(z,t) and 9{t) denote 
the background amplitude (which can be approximated 
in the framework of the Thomas- Fermi (TF) approxima- 
tion) and phase respectively, while the unknown complex 
function v{z,t) will represent a dark soliton. Assum- 
ing that the condensate dynamics involves a fast relax- 
ation scale to the ground state, and that the dark soliton 
evolves on top of this ground state, we obtain the fol- 
lowing perturbed nonlinear Schrodinger equation for the 
dark soliton wave function [21) : 



idtv + -d^v 



^\^-l)v = Piv), 



(10) 



where we have used the rescalings t ^>- fit and z -^ y/Jiz- 
The total perturbation P{v) in Eq. ( 10 ) has the form: 



Piv) 



1 
2/i 



2(1- 



\')V{z)v+^dzV- 



2fi^{z)dtv 
(11 



zo{t) 



ainif{t')dt', 



(13) 



is the soliton center. The evolution of the parameter (p 
is described by the following equation [63H65] : 



dip} 1 

dt 2 cos^ if sin (p 



Re 



+ CK3 



P{v)dtv''dz 



(14) 



The integral in Eq. ( 14 ) involves three terms [cf. 
Eq. (Ill]: the first two terms (accounting for the hamil- 



tonian part of the perturbation) can readily be evalu- 
ated upon expanding the potential V{z) in Taylor series 
around the soliton's center, zq; on the other hand, the 
third term (accounting for the dissipative part of the per- 
turbation) can be evaluated by approximating the func- 
tion 7(2) by its Taylor expansion around the trap center 
(z = 0). This expansion reads: 



7(z) « 70 



-72-2^ +^iZ^ +76^'^, 



where the constant coefficients are given by: 

7o = 2a/(c2-fd), 

j2 = 2a{3c'-d)/{c' + df, 

74 = 2a{5c^ - lOc^d + d^)/{c^ + df , 

76 = 2a{lc^ - 35c^d + 21c^d^ - d^)/{c^ 



(15) 



(16) 



-dy. 



we 



Using, as an example, the values of a, c and d corre- 
sponding to T = 150 nK (cf. Table |!|, in Fig. [11 
compare the fit to 7(2) given by Eq. (l8|) with the approx- 
imate expansion of Eq. ( [l5| . More specifically, we focus 
on the spatial interval of "i?/2, R/2] (where R is the TF 
radius of the cloud) - see inset of Fig. 11 this interval 



11 



is particularly relevant for our analytical and numerical 
considerations (see below), as we are mainly concerned 
with solitons moving close to the trap center. As seen in 
the figure, in this regime, the approximate expression is 
almost identical to the more accurate one, thus justifying 
the degree of approximation used in Eq. (15 1. 

To this end, using the approximate expression of 



Eq. (15), calculation of the integral in Eq. (14 1 leads to 



the following result: 



dLp 
Itt 



1 dV 2 

- cos ip — h 7o ^ M cos f sm ip 

2 dz 6 

{n^ - 6) 
■ 72^ -iita.n(p 
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74 
76 



7r2(7^2_go) ^^^ 2 



360 

TT* (3l7r2 - 294 

2016 



/i tan ipsec ip 



■/itaniy9sec cp. 



(17) 



Next, combining Eqs. (17) and (13), we obtain the fol 



lowing effective equation of motion for the dark soliton 
center. 




1- 



dzQ 
~dt 



76 



360 

7r'*(3l7r2-294) 
2016 



dzp 
dt 



1- 



dzQ 



H -2 



M^18) 



Notice that a variant of Eq. (18), corresponding to 72 — 
74 = 76 = 0, was presented in Ref. [21, . where the func- 
tion 7(2:) was approximated by the constant value 70 
(which is close to the mean value of 7(2:) in the interval 
[—R/2, R/2] - see Fig. 11 ). Nevertheless, here we are go- 



ing to analyze the more general problem and investigate 
the dissipative dynamics of solitons taking into regard 
the effect of the spatially-dependent profile of 7(2:). 

Here we should note that there appears to be a singu- 
larity in the solutions of Ec^. (llS]), corresponding to ve- 
locity values dzo/dt = 1, for which Eq. (18) becomes in- 
valid. This is also consistent with the analysis of Ref. [H7] , 
where it is shown that formal perturbation theory fails 
for extremely shallow dark solitons (with phase angles 
ip w 7r/2). In any case, in the physically relevant scenar- 
ios that we consider in our simulations below, solitons 
are observed to decay at the rims of the condensate at 
times smaller than the one needed for the soliton velocity 
to become dzo/dt = 1. 



B. The equation of motion for the soliton center 



The equation of motion ( 18 ) is evidently nonlinear 



Nevertheless, assuming that the dark soliton is close to 
a black one (i.e., (p is sufficiently small), Eq. (17) can be 
reduced to the following linearized form, 



'zo 



dt^ 



7oAi- 



dzQ 
dF 



n 



Zq ^ 0, 



(19) 



where we have introduced the variable 



70 



r°+- 18 

7r2 (7^2 _ 60) 
360" 



6) 



72 



-74 



TT^ (3l7r2 - 294) 
2016 



76. (20) 



Equation (19) is similar to the equation of motion de- 
rived in Ref. |23j by means of a kinetic-equation ap- 
proach. In the limiting case of zero temperature (i.e., 
for To = 0), Eq. (19) recovers the well-known result 

that a dark soliton 



(see, e.g., Refs. [15i [161 Hi 
oscillates with constant amplitude and frequency il/^/2 
in the harmonic trap V{z) — {l/2)il'^z'^. On the other 
hand, at finite temperatures (i.e., for 70 7^ 0), the lin- 
earized equation of motion ( |19[ ) additionally incorporates 
an anti-damping term, ex —dzg/dt [with a coefficient that 
takes into account the spatial dependence of 7(2)], which 
describes the expulsion of the dark soliton due to the 
interaction with the thermal cloud. 



It is clear that the nature of the solutions of Eq. ( 19 ) 



depend on whether the roots of the auxiliary equation 
s^ — (2/3)70^5 -I- (f7/\/2)2 = are real or complex. The 
roots are given as 



S1.2 = g7oM =^ ( 7s ) ^' ^ 



70^ 

7cr 



1, 



(21) 



where 7cr = (3/a*)(^/v2), and the discriminant A de- 
termines the type of the motion: 

In the super- critical case of strong anti-damping with 
A > 0, i.e., for high temperatures such that 70 > 7cr, 
the soliton trajectory is given by 



zo{t) 



1 



-[(s2Zo(0)-io(0))exp(sii) 
+ (io(0)-sizo(0))exp(s20]- 



(22) 



In the critical case with A 
trajectory is given as 



0, i.e.. To ~ 7cr, the soliton 



zo(t) = [zo(0) + (io(0) - ^^of^zo{0))t] expi^jof^t). (23) 

Finally, in the sub- critical case of weak anti-damping with 
A < 0, i.e., for sufficiently low temperatures such that 
To < 7cr, the soliton trajectory is given by 



Zo{t) =exp(-7o^i)[zo(0)cos(wosci) 

+ ^^03^(^0(0) - 7oM^o(0)) sm{ujosct)], 



(24) 
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FIG. 12. (Color online) Spatiotemporal evolution of the den- 
sity of a BEC confined in an harmonic trap with Q = 0.05 
and fi — 1, with a dark (black) soliton initally placed at 
zo = 2 for the super-critical case with jo/jcr = 1.89 (top 
panel; a — 84.44, c = 21.67 and d — 93.97), critical case 
with 7o/7cr ~ 1 (middle panel; a — 43.67, c = 22.19 and 
d = 59.16), and sub-critical case with 7o/7cr = 0.09 (bottom 
panel; a = 4.53, c = 23.89 and d = 19.72). Here -jcr = 0.106. 
The white dashed lines across the soliton trajectories corre- 
spond to the analytical prediction of the linearized equation 
of motion (191, while the black dashed lines are obtained by 
solving rmmorically the nonlinear equation of motion |T8| . 



where 



Wo 



n 




(25) 



is the soliton oscillation frequency. 

The above simple analysis shows that in the case of 
relatively high-temperatures (with 70 > 7cr), the dark 
soHton will not 'survive' long enough to oscillate in the 
trap, a result being in qualitative agreement with ex- 
perimental observations of dark matter-wave solitons in 
elongated Bose gases [T]. On the other hand, when the 
temperature is relatively low (i.e., 70 < "fcr), the dark 



soliton performs oscillations with an increasing amplitude 
and period - recall tha t the oscillation frequency is down- 
shifted as per Eq. (25) with respect to its value fl/\/2, at 
zero-temperature. These analytical results are in qual- 
itative agreement with the numerical ones obtained re- 
cently in the framework of the Zaremba-Nikuni-Grifhn 
model [milD]. 



C. Non-linear vs. linearized equation of motion 



To assess the above model, we now compare analytical 
[Eqs. ([l8l) - ([20|] and numerical [Eqs. ^ with 



soli- 
ton trajectories, in each case making use of the approxi- 
mate form for 7(2) employed in the analytical approach. 
First, we consider Eq. (|9]), with 7(2) given in Eq. ( [l5| , 
with an initial condition corresponding to a dark (black) 
soliton, initially placed off-centre at Zo{0) = 2 (with zero 
initial velocity, dzo{0)/dt = 0), on top of a TF cloud, 
with fi — 1, confined in a trap of strength fl = 0.05. 
The resulting soliton trajectories, found by integrating 
Eq. ([9]) by means of the split-step Fourier method, are 
shown in Fig. 12 for the super-critical (7/7^ — 1.89, top 
plot), critical (7 = jcr = 0.106, middle) and sub-critical 
filler = 0.09, bottom) cases. The black dashed lines 
in each plot depict the numerically obtained solutions of 
Eq. (18), while the white dashed lines show the corre- 



sponding analytical solutions of Eq. (19) 



Generally, it is observed that the results based on our 
analytical approximations, i.e., the solutions of Eqs. (18) 
and ( 19 ) - which were derived by employing the approx- 
imate form of 7(2) [cf. Eq. (11)] - are in good agree- 



ment with the numerical results obtained by the DGPE 
with the exact form of 7(2) [cf. Eq. (15)]. Neverthe- 
less, the solutions of the nonlinear equation of motion 



( 18 ) are more accurate than the analytical solutions of 



Eq. ( 19 ) in capturing the soliton trajectories obtained by 



the DGPE. This behavior is more pronounced for longer 
times, where the soliton either decays (top and middle 
panels of Fig. 12) or performs large amplitude oscilla- 
tions (bottom panel of Fig. 12): in fact, the solutions 



of the nonlinear equation of motion are able to correctly 
predict the decay time of the solitons [which is underes- 
timated by the solutions of Eq. ( 19 )] in the super-critical 



and critical cases of strong dissipation, or follow quite 
accurately the DGPE trajectory in the sub-critical case 
of weak dissipation; notice that, in the latter case, the 
analytical solution of Eq. ( 19 ) underestimates (overesti- 
mates) the frequency (amp. 
times. 



itude) of oscillation for longer 



We now proceed with a systematic comparison between 
analytical approximations, focussing on the more accu- 
rate nonlinear equation of motion (18), and numerical 
(DGPE and SGPE) resuhs. 
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FIG. 13. (Color online) Soliton trajectories as found by 
the numerical integration of the DGPE [solid (black) line], 
the SGPE [dash-dotted (red) line; single run with a decay 
time from the mean bin], and the nonlinear ordinary dif- 
ferential equation (ODE), Eq. ( |18[ ) [dashed (green) line], for 
T = 150 uK. Here a = 0.0219, c = 26.63, and d = 3.128. 



D. Comparison between numerical (SGPE and 
DGPE) results and analytical approximations 



To ensure long soliton lifetimes for this comparison, we 
focus on one of the lower temperatures within our study, 
T = 150 nK, which nevertheless still corresponds to 
ksT = 0.8/i. We perform a direct numerical integration 
of the DGPE, Eq. ^ [with j{z) given by Eq. ^] and 
compare this to respective results obtained via the ana- 
lytic ( 18 ) and SGPE models. The initial condition takes 



the form of a dark soliton, initially placed at the trap cen- 
tre Zo{0) = with initial velocity dzQ{0)/dt = 0.25 (the 
other parameters also remain the same with /i = l.58huj± 
and the dimensionless trap strength ft — 4 x 10^"^). 

In Fig. |13[ we show the soliton trajectories found via 
the DGPE and the SGPE, as well the solutions of the 



nonlinear equation of motion ( 18 ) stemming from our an- 



alytical approximations. It is clear that not only the solu- 
tion of the DGPE captures quite accurately the one of the 
SGPE (similarly to the behaviour found in Ref. [H]), but 
also the result of the analytical approximation is in very 
good agreement with the numerical results of the DGPE 
and SGPE. The SGPE trajectory shown is from a single 
run with a decay time close to the ensemble average; for 
these parameters, such trajectories were found to display 
dynamics close to the ensemble average in Ref. fIT\ . 

A further comparison between the analytical model 
above and the SGPE/DGPE can be seen in Fig. 10 where 
the predicted average decay times are plotted. In order 
to make this comparison, the decay times in the analyt- 
ical model must be extracted based on the visibility of 
the soliton over the background thermal fluctuations, as 
we discuss in the following Section. 



V. SOLITON VISUALIZATION 

A. Single shot fluctuation issues 

A quantity of relevance to experiments is the so-called 
visibility of the soliton. This is defined as ^Tlj : 



V 



(26) 



where rimax is the maximum BEG density and nmin is the 
minimum one, as set by the presence of the dark soliton. 
This parameter is a measure of how clearly a soliton can 
be seen in experiments. The SGPE accounts for both dis- 
sipation and fluctuations and can therefore be expected 
to produce experimentally-realistic visibility predictions; 
two examples of the visibility of solitons produced by the 
SGPE are shown in Fig. 14 (top and middle plots), show- 
ing an oscillatory decrease up to the point (denoted by 
the vertical dashed lines) when the soliton can no longer 
be distinguished from the background density fluctua- 
tions. 

In the displayed example trajectories, the solitons are 
lost to the background when the visibility decreases be- 
low around 50%. For comparison, we note that in the 
Hannover experiment [1], they reported a contrast in 
the range 20% — 40% throughout their measurements, 
and a visibilty of 50% in our present work corresponds 
to a contrast of nmin/nmax ^ 33%. The corresponding 
DGPE prediction is also shown (red curves) revealing 
good agreement in the region where the soliton can ac- 
tually be monitored over fluctuations in the stochastic 
cases (top and middle plots) . Clearly the results become 
meaningless beyond that time, as the soliton is then lost 
within that stochastic realization; betind this time the 
DGPE soliton signal remains visible, but this is due to 
incorrectly neglecting fluctuations. 

The periodic behaviour in the DGPE and SGPE vis- 
ibility arises because the soliton depth oscillates as the 
soliton traverses the harmonic trap. The points at which 
V = 1 correspond to the turning points of the soliton mo- 
tion, when the soliton depth is equal to the background 
density (or equivalently n,„in = 0). Perhaps a better 
measure is to look at the evolution at a specific point in 
the trap, e.g. the trap centre, which corresponds to the 
minimum visibility during the dissipative soliton motion. 
This is shown in Fig. 14 (bottom) and reveals quite good 
agreement between DGPE and SGPE, up to around the 
time that, on average, the SGPE solitons were visible. 

We now want to consider the corresponding analytical 
prediction based on our model of Section |IV[ As when 
comparing the numerical DGPE results to those of the 
SGPE (see Fig. 10 of Section HID), we again define the 



condition for soliton decay based upon the level of back- 
ground thermal fluctuations obtained from the SGPE. 
Then, calculating an expression for the soliton visibil- 
ity analytically, we are able to predict analytically the 
lifetime of a soliton within a finite temperature gas, ac- 
counting for both dissipation and background fluctua- 
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tions. Since ri,„in ^ fj, — /icos^(/9 (recall that the soli- 
ton depth is y/Ji cos (fi), the visibility can be expressed in 
terms of the soliton's phase angle as 



V = 



cos^ ip 



1 



• 2 

sm ip 



(27) 



Thus, < V < 1, with the limiting values V = and 
V = 1 corresponding, respectively, to a shallow soli- 
ton with if — >■ 7r/2 and a stationary kink with </?—>■ 0. 
Apparently, since ip ~ fit), the visibility is generally 
a function of time, but its analytical form can be de- 
termined via the time-dependence of (p{t), which can 
be derived numerically by means of Eq. (171 (see See- 
Nevertheless, it should be noticed that a 



IVAl 



tion 

simple analytical expression for the visibility can also 
be obtained in the case of sufficiently deep solitons 
(cos ip Ki 1 and sin ip k^ (jj) oscillating in a small re- 
gion around the trap center [i.e., for V{z = 0) = 0]: 
in this case, Eq. (17) becomes dp/dt = jof, and leads 
to the result p{t) — ipo{t) = iy9(0) exp(7o/ii)] (here, i^sq 
is the initial value of the phase angle) and, accordingly, 
VU=o = cos2 'Po{t)/[l + sin^ p^{t)]. 

Following the above arguments, we may estimate the 
soliton lifetime and the relevant results of the semi- 
analytical approximation, Eq. (27), are shown in Fig. 10 
of Section UTrni 

We can also calculate the soliton visibility versus 
time analytically, an example of which is shown in Fig- 
ure 14 [c) , alongside the numerical DGPE and average 
SGPE results at the same temperature. The SGPE re- 
sults in this case are an average over the visibility from 
few hundred runs, which smears out the oscillatory be- 
haviour, but yields a bulk behaviour consistent with that 
of the single realizations above. The agreement is very 
good between all approaches during the period that the 
soliton is visible over the background noise, beyond which 
time (vertical dashed lines), the SGPE visibility (black 
circles) plateaus. Here we have chosen to compare the 
analytical data to the minima of both the DGPE and 
average SGPE oscillations. A departure of the numeri- 
cal DGPE results from the SGPE data occurs close to 
the average soliton decay time, at which point many of 
the SGPE solitons have decayed, leaving a visibility read- 
ing in many single runs which corresponds to the back- 
ground noise signal. The DGPE soliton signal instead 
persists for much lower visibility values as the physical 
effects of fluctuations are not included, as also evident 
from the comparison of stochastic and dissipative 'car- 
pet plots' of Fig. [T] The analytical results are based on 
the Taylor expansion for "f{z) given by Eq. (15), which 



for large z is somewhat smaller than the full form of 7(2:) 
used in the numerical DGPE and SGPE simulations (see 
Fig. 11 ). This difference becomes apparent for times close 



to and after the average decay time (vertical dashed line) , 
for which the numerical and analytical results deviate 
slightly for this reason. 

Nonetheless, based on the level of background noise 
within the SGPE, and the analytical formula for the vis- 



SGPE (single run) 
DGPE 




FIG. 14. (Color online) Numerical results for the visibility as 
a function of time for two example SGPE single runs (noisy 
turquoise, top plot; noisy brown, middle plot) vs. the DGPE 
(smooth red curve in top and middle plots); the insets show 
the corresponding soliton trajectories. The lower plot shows 
the analytical result [green squares] versus the average visi- 
bility from the SGPE (black circles) and the visibility from 
the DGPE (red diamonds; the latter two data sets correspond 
to the local minima of oscillations, which can be seen in the 
upper plots. At this temperature {T = 175nK), the analyt- 
ical prediction is very good up to the point that the soliton 
can be tracked, though uses the approximate form for 7(2) of 
Eq. ^ while the SGPE and DGPE use Eq.^. The stochas- 
tic results are examples of solitons with decay times close to 
the ensemble mean, indicated by the vertical dashed line in 
all plots. 



ibility, it should be possible to predict realistic lifetimes 
for solitons within experiments based on this approach. 



B. Comments on related work 

The study of dark solitons using stochastic and clas- 
sical field methods has received significant attention re- 
cently pn 13^17^75] . In particular, the SGPE was ap- 
plied initially in parallel, independent works [211 EH] , and 
recently also in [75]. The analysis presented in TT con- 
sidered dark solitons as relics due to a quench of the sys- 
tem parameters within a homogeneous, periodic ID Bose 
gas. Notably, the distribution reported in ^3\ describing 
the number of solitons with time (Fig. 3 of Ref. [TSj) has 
a form which is qualitatively similar to the decay time 
distributions of Refs. ^2ll |55] and in the present work. 

More recently, Wright and Bradley [76] performed a 
study which is more closely related to Ref. [21] and our 
current analysis, based on the stochastic projected Gross- 
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Pitaevskii equation (SPGPE) [HllZZj- (This method is a 
variant of the SGPE ,38^ which features a projector into 
low energy modes.) They found that the stochastic sim- 
ulations yielded average velocities which were lower than 
those found within dissipative GPE simulations, imply- 
ing that the solitons have longer lifetimes, on average, 
within the stochastic simulations, i.e., that the noise pro- 
longs the solitons' existence. This is in agreement with 
the long tails in the decay time distributions of both 
Refs. [211 [55] and the current work. Furthermore, the 
analytical findings of Ref. [75], and in particular Eq.(18) 
of Ref. [75], arise as a special case of our previously re- 
ported results [U (see sentence preceeding Eq. (6) in 
[21]) and by extension also of this work, which addition- 
ally treats a spatially-dependent dissipative term. This 
can be seen from Eq. (17) above by setting V — (for 
a homogeneous system) and replacing 7(2) —>■ 70 (i.e. 
72 = 74 = 76 = for a spatially-constant dissipation), 
and multiplying through by cos{(j)) [78] . 

Beyond this, it is not straightforward to give a more 
direct comparison to their findings for several reasons: (i) 
we consider a harmonically trapped system, while they 
focussed on soliton propagation within a homogeneous 
sample; (ii) they measure the increase in velocity as the 
soliton decays, which is straightforward without a trap, 
while in our case the soliton velocity is constantly chang- 
ing as it oscillates; an upshot of this is that the soliton 
depth does not change monotonically as it decays, but 
has an additional oscillatory component, complicating 
the velocity analysis; (iii) we analyze the ensemble of soli- 
tons via their decay times, and therefore effectively sam- 
ple members of the SGPE ensemble at different times, 
while they consider equal time measurements of the ve- 
locity via an ensemble measurement at a particular time. 



VI. CONCLUSIONS 

We have characterized the dynamics of dark solitons 
propagating within a partially condensed, harmonically 
trapped Bose gas in the presence of phase fluctuations. 
Our analysis was based on a stochastic Gross-Pitaevskii 
equation, which reduces to a dissipative Gross-Pitaevskii 
equation upon neglect of the additive noise term. 

Stochastic simulations allowed us to perform a statis- 
tical analysis on the soliton decay times. Our results 
showed that to study soliton decay, information should 
be extracted from single soliton realizations prior to per- 
forming such an analysis, in agreement with related ex- 
perimental findings [5D]. On doing so, we found dark 
soliton lifetimes to be approximately lognormally dis- 
tributed, which implies that some solitons within a num- 
ber of realizations may have very long lifetimes relative 
to the ensemble mean, an effect already observed in ex- 



periments [5]. We found the standard deviation and 
skewness of these distributions to increase monotonically, 
and approximately linearly with temperature. Extract- 
ing expectation values from the decay time distribution 
obtained at each temperature, showed the purely dissipa- 
tive results matched these stochastic expectation values 
decay times well, once the effects of background fluctua- 
tions were taken into account. 

Considering the interplay between noise in the initial 
conditions (used e.g. in some simpler approximate mod- 
els) and the dynamical noise at each temporal step of 
the Stochastic Gross-Pitaevskii equation - which loosely 
corresponds to a stochastic random kick to the soliton 
position -, we found the dynamical noise to play an im- 
portant role in determining the final decay time. 

We also presented results for the experimentally rele- 
vant visibility of the soliton and found the soliton decay 
to be related to the visibility in our numerical simulations 
reaching a plateau (at which point our soliton tracking 
algorithm breaks down) . The value of the plateau is set 
by the strength of the background fluctuations, which is 
a direct measure of temperature in these systems. 

In the purely dissipative case, we derived analytical 
expressions describing the dynamics of the soliton den- 
sity notch, with good agreement found between these 
and numerical solution to the dissipative Gross-Pitaevskii 
equation, generalizing our earlier work [21] to the case 
of a spatially dependent damping, as obtained ab initio 
within the stochastic Gross-Pitaevskii formalism. The 
average soliton decay times were found to scale as T""*; 
this has been previously obtained for homogeneous sys- 
tems at low temperatures ksT <C ^J-, but our numerics 
indicates that (at least within the classical field approxi- 
mation for the low-lying modes of the system) this can be 
extended to the trapped case, and to the regime for which 
kgT < /i, even in the presence of phase fiuctuations. 

Observing the dynamics of dark solitons within phase- 
fluctuating condensates offers an intriguing opportunity 
to observe a macroscopic quantum object undergoing a 
Brownian-like motion. We hope that the fluctuating as- 
pects of this behavior, and especially the dependence on 
temperature, may be further analyzed in future experi- 
ments on ultracold gases in the near future. 
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